package ij3d.Jama;

import java.io.Serializable;

/**
 * Cholesky Decomposition.
 *<P>
 *For a symmetric, positive definite matrix A, the Cholesky decomposition
 *is an lower triangular matrix L so that A = L*L'.
 *<P>
 *If the matrix is not symmetric or positive definite, the constructor
 *returns a partial decomposition and sets an internal flag that may
 *be queried by the isSPD() method.
 *
 * @author     thomas
 * @created    18 décembre 2007
 */

public class CholeskyDecomposition implements java.io.Serializable {

	/* ------------------------
	 *Class variables
	 * ------------------------  */
	/**
	 * Array for internal storage of decomposition.
	 *
	 * @serial    internal array storage.
	 */
	private double[][] L;

	/**
	 * Row and column dimension (square matrix).
	 *
	 * @serial    matrix dimension.
	 */
	private int n;

	/**
	 * Symmetric and positive definite flag.
	 *
	 * @serial    is symmetric and positive definite flag.
	 */
	private boolean isspd;


	/* ------------------------
	 *Constructor
	 * ------------------------  */
	/**
	 * Cholesky algorithm for symmetric and positive definite matrix.
	 *
	 * @param  Arg  Description of the Parameter
	 */

	public CholeskyDecomposition(Matrix Arg) {
		// Initialize.
		double[][] A = Arg.getArray();
		n = Arg.getRowDimension();
		L = new double[n][n];
		isspd = (Arg.getColumnDimension() == n);
		// Main loop.
		for (int j = 0; j < n; j++) {
			double[] Lrowj = L[j];
			double d = 0.0;
			for (int k = 0; k < j; k++) {
				double[] Lrowk = L[k];
				double s = 0.0;
				for (int i = 0; i < k; i++) {
					s += Lrowk[i] * Lrowj[i];
				}
				Lrowj[k] = s = (A[j][k] - s) / L[k][k];
				d = d + s * s;
				isspd = isspd & (A[k][j] == A[j][k]);
			}
			d = A[j][j] - d;
			isspd = isspd & (d > 0.0);
			L[j][j] = Math.sqrt(Math.max(d, 0.0));
			for (int k = j + 1; k < n; k++) {
				L[j][k] = 0.0;
			}
		}
	}


	/* ------------------------
	 *Temporary, experimental code.
	 * ------------------------ *\
	 *\** Right Triangular Cholesky Decomposition.
	 *<P>
	 *For a symmetric, positive definite matrix A, the Right Cholesky
	 *decomposition is an upper triangular matrix R so that A = R'*R.
	 *This constructor computes R with the Fortran inspired column oriented
	 *algorithm used in LINPACK and MATLAB.  In Java, we suspect a row oriented,
	 *lower triangular decomposition is faster.  We have temporarily included
	 *this constructor here until timing experiments confirm this suspicion.
	 *\
	 *\** Array for internal storage of right triangular decomposition. **\
	 *private transient double[][] R;
	 *\** Cholesky algorithm for symmetric and positive definite matrix.
	 *@param  A           Square, symmetric matrix.
	 *@param  rightflag   Actual value ignored.
	 *@return             Structure to access R and isspd flag.
	 *\
	 *public CholeskyDecomposition (Matrix Arg, int rightflag) {
	 */ 
	 /*
	 Initialize.
	 *double[][] A = Arg.getArray();
	 *n = Arg.getColumnDimension();
	 *R = new double[n][n];
	 *isspd = (Arg.getColumnDimension() == n);
	 */ 
	 /*Main loop.
	 *for (int j = 0; j < n; j++) {
	 *double d = 0.0;
	 *for (int k = 0; k < j; k++) {
	 *double s = A[k][j];
	 *for (int i = 0; i < k; i++) {
	 *s = s - R[i][k]*R[i][j];
	 *}
	 *R[k][j] = s = s/R[k][k];
	 *d = d + s*s;
	 *isspd = isspd & (A[k][j] == A[j][k]);
	 *}
	 *d = A[j][j] - d;
	 *isspd = isspd & (d > 0.0);
	 *R[j][j] = Math.sqrt(Math.max(d,0.0));
	 *for (int k = j+1; k < n; k++) {
	 *R[k][j] = 0.0;
	 *}
	 *}
	 *}
	 *\** Return upper triangular factor.
	 *@return     R
	 *\
	 *public Matrix getR () {
	 *return new Matrix(R,n,n);
	 *}
	 *\* ------------------------
	 *End of temporary code.
	 * ------------------------  */
	/* ------------------------
	 *Public Methods
	 * ------------------------  */
	/**
	 * Is the matrix symmetric and positive definite?
	 *
	 * @return    true if A is symmetric and positive definite.
	 */

	public boolean isSPD() {
		return isspd;
	}


	/**
	 * Return triangular factor.
	 *
	 * @return    L
	 */

	public Matrix getL() {
		return new Matrix(L, n, n);
	}


	/**
	 * Solve A*X = B
	 *
	 * @param  B                             A Matrix with as many rows as A and any number of columns.
	 * @return                               X so that L*L'*X = B
	 */

	public Matrix solve(Matrix B) {
		if (B.getRowDimension() != n) {
			throw new IllegalArgumentException("Matrix row dimensions must agree.");
		}
		if (!isspd) {
			throw new RuntimeException("Matrix is not symmetric positive definite.");
		}

		// Copy right hand side.
		double[][] X = B.getArrayCopy();
		int nx = B.getColumnDimension();

		// Solve L*Y = B;
		for (int k = 0; k < n; k++) {
			for (int i = k + 1; i < n; i++) {
				for (int j = 0; j < nx; j++) {
					X[i][j] -= X[k][j] * L[i][k];
				}
			}
			for (int j = 0; j < nx; j++) {
				X[k][j] /= L[k][k];
			}
		}

		// Solve L'*X = Y;
		for (int k = n - 1; k >= 0; k--) {
			for (int j = 0; j < nx; j++) {
				X[k][j] /= L[k][k];
			}
			for (int i = 0; i < k; i++) {
				for (int j = 0; j < nx; j++) {
					X[i][j] -= X[k][j] * L[k][i];
				}
			}
		}
		return new Matrix(X, n, nx);
	}
}

